from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
import numpy as np
import matplotlib.cm as cm
import matplotlib.colors as cls
from pandas import Series, DataFrame
from mpl_toolkits.axes_grid1.inset_locator import inset_axes, zoomed_inset_axes, mark_inset
import matplotlib.ticker as ticker
from collections import defaultdict, namedtuple
from pandas import ExcelWriter
import xlsxwriter
import pandas as pd
import calendar
import pickle
import datetime
from scipy import stats
import itertools
import statistics
# set seaborn plot style
import seaborn as sns
sns.set()
sns.set_context("paper")
sns.set_style("whitegrid")
# set default font
plt.rcParams["font.family"] = "Times New Roman"
def get_ecos():
m = Basemap(projection='cyl')
shpf = m.readshapefile('C:/Users/alpha/Documents/MODIS/MCD64A1/newreprojected/newreprojected','ecos',linewidth=0.1)
return m.ecos_info, m.ecos
m = Basemap(projection='cyl')
shpf = m.readshapefile('C:/Users/alpha/Documents/MODIS/MCD64A1/newreprojected/newreprojected','ecos',linewidth=0.1)
ecos_info, ecos = m.ecos_info, m.ecos
%qtconsole
ecos_info[0].keys()
months = [calendar.month_name[i +1][:3] for i in range(12)]
years = [str(y) for y in range(2001,2018)]
ecoids = [e['ECO_ID'] for e in ecos_info]
bionums = [e['BIOME_NUM'] for e in ecos_info]
bionames = [e['BIOME_NAME'] for e in ecos_info]
ids = [e['ID'] for e in ecos_info]
climes = [e['GRIDCODE'] for e in ecos_info]
continents = [e['CONTINENT'] for e in ecos_info]
countries = [e['FIPS_CNTRY'] for e in ecos_info]
areas= [e['Area_polig'] for e in ecos_info]
for i ,id in enumerate(ecoids):
if id == 0:
bionums[i] = 0.0
result = {}
for i in range(1,18):
suf = '20{0:02d}'.format(i)
in_file= 'newresult' + suf
result[suf] = DataFrame(np.load(in_file), columns =months)
for i in range(1,18):
suf = '20{0:02d}'.format(i)
result[suf] = result[suf].assign(eco_id=ecoids,biome_num =bionums,biome_name=bionames,idd=ids,clime=climes,continent=continents,country=countries, area = areas)
grouped = result['2001'][ ['eco_id','continent','clime','country','idd' ,'area'] ].groupby(['eco_id','continent','clime','country','idd' ,'area'])
areas_row = []
for a, _ in grouped:
areas_row.append(a[-1])
attrs_table = {}
for yr, df in result.items():
attrs_table[yr] = df[ months + ['eco_id','continent','clime','country','idd','biome_num' ,'area']].groupby(['eco_id','continent','clime','country','idd','biome_num','area']).sum()
at_table = pd.concat(attrs_table)
writer = ExcelWriter('newshape.xlsx',engine='xlsxwriter')
at_table.to_excel(writer,sheet_name='Burned Pixels',engine='xlsxwriter')
writer.save()
plt.rcParams['mathtext.fontset'] = 'custom'
plt.rcParams['mathtext.rm'] = "Times New Roman"
plt.rcParams['mathtext.it'] = "Times New Roman"
plt.rcParams['mathtext.bf'] = "Times New Roman"
g = sns.regplot(np.arange(12),attrs_table['2002'].sum()/1e6, fit_reg= False)
ylabels = g.get_yticklabels()
g.set_xticklabels(months[::2], fontsize=15)
g.set_xticks(np.arange(0,12,2))
g.set_ylabel('Burned Pixels $(10^6)$', fontsize=15)
plt.setp(ylabels,fontsize=15)
#g.set(xticklabels=months[::2],xticks=np.arange(0,12,2), ylabel='Burned Pixels $(10^6)$', ylabelsize=15)
plt.show()
plt.plot(np.arange(12),attrs_table['2002'].sum()/1e6)
g =plt.gca()
ylabels = g.get_yticklabels()
g.set_xticklabels(months[::2], fontsize=15)
g.set_xticks(np.arange(0,12,2))
g.set_ylabel('Burned Pixels $(10^6)$', fontsize=15)
plt.setp(ylabels,fontsize=15)
#g.set(xticklabels=months[::2],xticks=np.arange(0,12,2), ylabel='Burned Pixels $(10^6)$', ylabelsize=15)
plt.show()
plt.scatter(np.arange(12),attrs_table['2002'].sum()/1e6)
g =plt.gca()
ylabels = g.get_yticklabels()
g.set_xticklabels(months[::2], fontsize=15)
g.set_xticks(np.arange(0,12,2))
g.set_ylabel('Burned Pixels $(10^6)$', fontsize=15)
plt.setp(ylabels,fontsize=15)
#g.set(xticklabels=months[::2],xticks=np.arange(0,12,2), ylabel='Burned Pixels $(10^6)$', ylabelsize=15)
plt.show()
at_table.shape
at_table.loc['2001',0,'Antarctica']
at_table.sum()
fig, ax = plt.subplots()
ax.plot( np.arange(12),at_table.sum()/1e6)
ax.set_xticks(np.arange(12))
ax.set_xticklabels(months, fontsize=15)
ax.set_title('Intra-annual Variability',fontsize=15)
labels_x = ax.get_xticklabels()
labels_y = ax.get_yticklabels()
off =ax.get_yaxis().get_offset_text()
off.set_x(-.07)
off.set_fontsize(12)
plt.setp(labels_y,fontsize=15)
plt.ylabel(r"Burned Pixels $(10^6)$", fontsize = 15)
plt.tight_layout()
plt.show()
at_table.index.names = ['year', 'eco_id', 'continent', 'clime', 'country', 'idd', 'biome_num', 'area']
at_table.groupby('year').sum().sum(axis=1)
fig, ax = plt.subplots()
labels = years[::2]
ax.plot(years, at_table.groupby('year').sum().sum(axis=1)/1e6)
ax.set_xticks(labels)
ax.set_xticklabels(labels, fontsize=15)
ax.set_title('Interannual Variability',fontsize=15)
labels_x = ax.get_xticklabels()
labels_y = ax.get_yticklabels()
off =ax.get_yaxis().get_offset_text()
off.set_x(-.07)
off.set_fontsize(12)
plt.setp(labels_y,fontsize=15)
plt.ylabel(r"Burned Pixels ($10^6$)", fontsize = 15)
plt.tight_layout()
plt.show()
at_table.groupby('year').sum().values.flatten()
ym = itertools.product(years,months)
ymlabels = [ m + ' ' + y[-2:] for y, m in ym]
fig, ax = plt.subplots(figsize=(16,12))
labels = ymlabels[::6]
ax.plot( np.arange(12*17),at_table.groupby('year').sum().values.flatten()/1e6)
ax.axhline(y=statistics.mean(at_table.groupby('year').sum().values.flatten()/1e6), color='r', linestyle='--', label='mean')
ax.set_xticks(np.arange(0,12*17,6))
ax.set_xticklabels(labels,fontsize=12, rotation=30)
labels_y = ax.get_yticklabels()
plt.setp(labels_y,fontsize=15)
plt.ylabel(r"Burned Pixels ($10^6$)", fontsize = 15)
ax.set_title('Intra-monthly Variability',fontsize=15)
plt.legend(fontsize=15)
#prop={'size': 6}
plt.show()
at_table.groupby('biome_num').sum()
at_table.loc['2001',0, 'Antarctica',61,'AY',2250 ]
areas = at_table.loc['2001'].index.get_level_values(6)
at_table.index.get_level_values(7)
141338/17
all_areas =at_table.index.get_level_values(7)
at_table_cp = at_table.copy()
#at_table_cp.areas = all_areas
at_table_cp.reset_index(inplace=True)
agg_dict = {month : 'sum' for month in months}
div_17 = lambda x: sum(x)/17
agg_dict['area'] = div_17
biomes_table = at_table_cp.groupby('biome_num').agg(agg_dict)
biomes_table[months].div(biomes_table['area'], axis =0)
biomes_table['area']
# for serialization:
tabelas = {'multindex': at_table, 'simpleindex':at_table_cp}
with open('tabelas','wb') as f:
f.write(pickle.dumps(tabelas))
biomes_intra = biomes_table[months].div(biomes_table['area'], axis =0)
biome_names_dict = dict( (num, name) for num, name in zip(bionums, bionames))
fig, ax = plt.subplots(5,3, figsize =(18,12) , sharey=True)
labels = months
ax = ax.ravel()
for i in range(15):
#slope, intercept, r_value, p_value, std_err = stats.linregress(np.arange(12), monbio[i*3 + j])
#ax[i][j].plot(months, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i].plot(months, biomes_intra.iloc[i], 'o')
ax[i].set_xticks(np.arange(12))
ax[i].set_xticklabels(labels,fontsize=15)
labels_y = ax[i].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i].set_title(biome_names_dict[i],fontsize=15)
#ax[i].legend(fontsize=15)
fig.tight_layout()
plt.show()
fig, ax = plt.subplots(5,3, figsize =(18,12))
labels = months
ax = ax.ravel()
for i in range(15):
#slope, intercept, r_value, p_value, std_err = stats.linregress(np.arange(12), monbio[i*3 + j])
#ax[i][j].plot(months, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i].plot(months, biomes_intra.iloc[i], 'o')
ax[i].set_xticks(np.arange(12))
ax[i].set_xticklabels(labels,fontsize=15)
labels_y = ax[i].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i].set_title(biome_names_dict[i],fontsize=15)
#ax[i].legend(fontsize=15)
fig.tight_layout()
plt.show()
at_table_cp[at_table_cp.year == '2001'].groupby('biome_num').sum().div(biomes_table['area'], axis=0)[months]
biomes_table['area'].size
## test 2001 may bionum = 5 assert 3.547547e-04
bp = at_table_cp[at_table_cp.year == '2001'].groupby('biome_num').sum().iloc[5]['May']
area_5 = biomes_table['area'][5]
assert bp/area_5 == 3.547547e-04, print((bp/area_5 - 3.547547e-04 )*100/3.547547e-04)
at_table_cp[at_table_cp.year == '2001'].groupby('biome_num').sum().div(biomes_table['area'], axis=0)[months].sum(axis=1)
biomes_inter = [at_table_cp[at_table_cp.year == yr].groupby('biome_num').sum().div(biomes_table['area'], axis=0)[months].sum(axis=1) for yr in years]
biomes_inter = DataFrame(biomes_inter).T
biomes_inter
fig, ax = plt.subplots(5,3, figsize =(18,12) , sharey=True)
labels = years[::2]
ax = ax.ravel()
for i in range(15):
slope, intercept, r_value, p_value, std_err = stats.linregress(range(17), biomes_inter.iloc[i])
ax[i].plot(years, intercept + slope*range(17), 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i].plot(years, biomes_inter.iloc[i], 'o', label='')
ax[i].set_xticks(labels)
ax[i].set_xticklabels(labels,fontsize=15)
labels_y = ax[i].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i].set_title(biome_names_dict[i],fontsize=15)
ax[i].legend(fontsize=15)
fig.tight_layout()
plt.show()
fig, ax = plt.subplots(5,3, figsize =(18,12) )
labels = years[::2]
ax = ax.ravel()
for i in range(15):
#slope, intercept, r_value, p_value, std_err = stats.linregress(range(17), biomes_inter.iloc[i])
#ax[i].plot(years, intercept + slope*range(17), 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i].plot(years, biomes_inter.iloc[i], label='')
ax[i].set_xticks(labels)
ax[i].set_xticklabels(labels,fontsize=15)
labels_y = ax[i].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i].set_title(biome_names_dict[i],fontsize=15)
ax[i].legend(fontsize=15)
fig.tight_layout()
plt.show()
biomes_inter
at_table_cp[at_table_cp.year == '2001'].groupby('biome_num').sum().div(biomes_table['area'], axis=0)[months].sum()
no_climas = len(set(climes));no_climas
climes_table = at_table_cp.groupby('clime').agg(agg_dict)
climes_table[months].div(climes_table['area'], axis =0)
climes_intra = climes_table[months].div(climes_table['area'], axis =0)
climanames = [e['Clima_Name'] for e in ecos_info]
climes_intra.index[0]
clima_names_dict = dict( (num, name) for num, name in zip(climes, climanames))
# 11 x 3, 31 climas
fig, ax = plt.subplots(nrows=11,ncols=3, figsize =(18,40) , sharey=True)
labels = months
ax = ax.ravel()
for i in range(31):
#slope, intercept, r_value, p_value, std_err = stats.linregress(np.arange(12), monbio[i*3 + j])
#ax[i][j].plot(months, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i].plot(months, climes_intra.iloc[i], 'o')
ax[i].set_xticks(np.arange(12))
ax[i].set_xticklabels(labels,fontsize=15)
labels_y = ax[i].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i].set_title(clima_names_dict[climes_intra.index[i]],fontsize=15)
#ax[i].legend(fontsize=15)
for i in range(31,33):
ax[i].set_visible(False)
fig.tight_layout()
plt.show()
fig, ax = plt.subplots(nrows=11,ncols=3, figsize =(18,40) )
labels = months
ax = ax.ravel()
for i in range(31):
#slope, intercept, r_value, p_value, std_err = stats.linregress(np.arange(12), monbio[i*3 + j])
#ax[i][j].plot(months, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i].plot(months, climes_intra.iloc[i], 'o')
ax[i].set_xticks(np.arange(12))
ax[i].set_xticklabels(labels,fontsize=15)
labels_y = ax[i].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i].set_title(clima_names_dict[climes_intra.index[i]],fontsize=15)
#ax[i].legend(fontsize=15)
for i in range(31,33):
ax[i].set_visible(False)
fig.tight_layout()
plt.show()
clima_names_dict.values()
at_table_cp[at_table_cp['year']=='2001'].groupby('clime').sum().div(climes_table['area'], axis=0)[months].sum(axis=1)
climas_inter = [at_table_cp[at_table_cp['year']==yr].groupby('clime').sum().div(climes_table['area'], axis=0)[months].sum(axis=1) for yr in years]
climas_inter = DataFrame(climas_inter).T
climas_inter
fig, ax = plt.subplots(nrows=11,ncols=3, figsize =(18,40) , sharey=True)
labels = years[::2]
ax = ax.ravel()
for i in range(31):
slope, intercept, r_value, p_value, std_err = stats.linregress(range(17), climas_inter.iloc[i])
ax[i].plot(years, intercept + slope*range(17), 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i].plot(years, climas_inter.iloc[i], 'o', label='')
ax[i].set_xticks(labels)
ax[i].set_xticklabels(labels,fontsize=15)
labels_y = ax[i].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i].set_title(clima_names_dict[climas_inter.index[i]],fontsize=15)
ax[i].legend(fontsize=15)
for i in range(31,33):
ax[i].set_visible(False)
fig.tight_layout()
plt.show()
fig, ax = plt.subplots(nrows=11,ncols=3, figsize =(18,40) )
labels = years[::2]
ax = ax.ravel()
for i in range(31):
slope, intercept, r_value, p_value, std_err = stats.linregress(range(17), climas_inter.iloc[i])
ax[i].plot(years, intercept + slope*range(17), 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i].plot(years, climas_inter.iloc[i], 'o', label='')
ax[i].set_xticks(labels)
ax[i].set_xticklabels(labels,fontsize=15)
labels_y = ax[i].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i].set_title(clima_names_dict[climas_inter.index[i]],fontsize=15)
ax[i].legend(fontsize=15)
for i in range(31,33):
ax[i].set_visible(False)
fig.tight_layout()
plt.show()
groupedby_year_clima = at_table_cp.groupby(['clime','year']).sum()
groupedby_year_clima = groupedby_year_clima[months].div(groupedby_year_clima.area, axis =0)
groupedby_year_clima =groupedby_year_clima.reset_index()
groupedby_year_clima[groupedby_year_clima['clime'] == 11][months].T.iloc[0]
labels = years[::2]
for climcode in clima_names_dict.keys():
df = groupedby_year_clima[groupedby_year_clima['clime'] == climcode][months].T
fig, ax = plt.subplots(4,3, figsize =(20,12) , sharey= True)
ax = ax.ravel()
fig.suptitle(clima_names_dict[climcode],fontsize=15)
for i in range(12):
data = df.iloc[i]
slope, intercept, r_value, p_value, std_err = stats.linregress(range(17), data)
ax[i].plot(years, intercept + slope*range(17), 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i].plot(years, data, 'o', label='_nolegend_')
ax[i].set_xticks(labels)
ax[i].set_xticklabels(labels,fontsize=15)
labels_y = ax[i].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i].set_title(calendar.month_name[i +1],fontsize=15)
ax[i].legend(fontsize=15)
fig.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
labels = years[::2]
for climcode in clima_names_dict.keys():
df = groupedby_year_clima[groupedby_year_clima['clime'] == climcode][months].T
fig, ax = plt.subplots(4,3, figsize =(20,12) )
ax = ax.ravel()
fig.suptitle(clima_names_dict[climcode],fontsize=15)
for i in range(12):
data = df.iloc[i]
#slope, intercept, r_value, p_value, std_err = stats.linregress(range(17), data)
#ax[i].plot(years, intercept + slope*range(17), 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i].plot(years, data, 'o', label='_nolegend_')
ax[i].set_xticks(labels)
ax[i].set_xticklabels(labels,fontsize=15)
labels_y = ax[i].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i].set_title(calendar.month_name[i +1],fontsize=15)
ax[i].legend(fontsize=15)
fig.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
# as linhas seguintes dão o mesmo resultado
at_table_cp.groupby(['year','clime']).sum()[months];
pd.pivot_table(at_table_cp, values=months, index=['year', 'clime'],aggfunc=np.sum)[months];
notincluded13 = pickle.loads(open('notincluded2013','rb').read())
df_notinc13 = DataFrame(notincluded13)
df_notinc13.to_csv('not_included13.csv', header=False)
notincluded17 = pickle.loads(open('notincluded2017','rb').read())
df_notinc17 = DataFrame(notincluded17)
df_notinc17.to_csv('not_included17.csv', header=False)